/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
/** @file cArrays.cpp
 *  @brief Implementation of cVector and cMatrix classes of complex numbers
 */

#include "cArrays.h"

using namespace std;
namespace smv_lib
{

//Defining ARRAYS_DEBUG turns on bounds checking and size checking
#define ARRAYS_DEBUG

/*==================================================
           cVector class definitions
 elements of cVector are complex<double> numbers
     Complex z  => z = x + iy   (i*i = -1)
==================================================*/
//--- constructors
cVector::cVector()                      
{
	data = new Complex[M_DEFAULT];
	assert(data != NULL);
	len = M_DEFAULT;
	Name = "cVECTOR";
	isNameGiven = false;
}

cVector::cVector(int n, Complex d) : len(n)
{
	data = new Complex[n];
	assert(data != NULL);
	for (int i = 0; i < n; i++)
		data[i] = d;
	Name = "cVECTOR";
	isNameGiven = false;
}

cVector::cVector(int n, string name) : len(n), Name(name)
{
	data = new Complex[n];
	assert(data != NULL);
	Complex d (0.0, 0.0);
	for (int i = 0; i < n; i++)
		data[i] = d;
	if (Name == "") {
		Name = "cVECTOR";
	    isNameGiven = false;
	}
	isNameGiven = true;
}

/**
** \brief Construct a cVector from a file
**
**    Format of the Vector file - it contains two parts: Header block and Data block.
**    HEADER BLOCK: 
**
**    (1) One or more optional COMMENT lines, where each comment line begins with '!'
**    
**    (2) Optional line with a NAME of the Vector. If given, it must be written as:
**    N= (string), 
**    where the string in the name must begin with a letter (a-w,A-W).
**    If the Name is not given, default name VECTOR is assigned.
**
**    (3) Number of elements 'n' and number of entries in a line 'l'
**    (for small vectors n and l are equal - all in one line).
**    Both entries 'n' and 'l' must be given!!
**    n= (int)  l= (int)
**           
**    The end of the Header block is given as 3 colons:
**    :::
**    DATA BLOCK:
**    numeric entris are given in a row 
**    (entries are separated by a blanko)
**
**    NOTE: For longer vectors (say, n > 15-20), multiline input is neccessary 
**    
**/
static char cAlphabet[] = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ";

cVector::cVector(const char* file)
{
    const int BSIZE = 120;
    FILE *from;
    int nele  = 0;  // total number of elements in a cVector
    int mult  = 0;  // number of elements in a single line

    if ((from = fopen(file,"r")) == 0) {
	cerr << "\nCannot open " << file << endl;
	exit(1);
    }
    char buf[BSIZE], *cp, *cp2;
    int efound = 0, mfound = 0, nfound = 0, colonsfound = 0;

    //--------------------------------------------- read the HEADER BLOCK
    while( fgets(buf, BSIZE, from) ) {
		//--- read the comment lines
      if( ( cp = strpbrk(buf,"!")) != NULL )         // comment line
        *cp = '\0';                                  // terminate string at comment
		//---- read the Name of a Vector 
      if(( cp = strpbrk(buf,"N")) != NULL )          
        if (( cp2 = strpbrk(cp, "=")) != NULL )
          if (( cp = strpbrk( cp2, cAlphabet)) != NULL ) {
			  Name = cp;
			  isNameGiven = true;
			  nfound++;
            }   
		  //--- read the lenth of a Vector
	  if(( cp = strpbrk(buf,"n")) != NULL )
        if (( cp2 = strpbrk(cp, "=")) != NULL )
          if (( cp = strpbrk( cp2, "0123456789")) != NULL ) {
              nele = atoi(cp);                       
              efound++;                              
            }
		  //--- read the number of entries in a line
      if(( cp = strpbrk(buf,"l")) != NULL )
        if (( cp2 = strpbrk(cp, "=")) != NULL )
          if (( cp = strpbrk( cp2, "0123456789")) != NULL ) {
              mult = atoi(cp);                       
              mfound++;                              
            }
		  //--- termination of the Header block
      if ( strstr(buf,":::") != NULL ) {             
          colonsfound++;
          break; //... out of "while" loop
        }
    }
    if ( !efound || ! mfound  || !colonsfound ) {
		cerr << "\nThe input file is not OK \n" << endl;
        exit(1);
    }
	if (!nfound) {               // no Vector name is given
		Name = "cVECTOR";         // assign the default name
		isNameGiven = false;
	}
	//------------------------------------------ make space for the cVector
	data = new Complex[nele];
	len  = nele;
	int i = 0;
	double num;
	// number of line to read 
	int No_Lines, resto;
	div_t Res;
	Res.quot = 1;  Res.rem = 0;
	if (nele > mult) 
		Res = div(nele, mult);
	No_Lines = Res.quot;        // No_Lines with full number of mult entries
	resto    = Res.rem;         // ONE line with resto entries

    //-------------------------------------------- read the DATA BLOCK    
	int kmax;
	for (int j = 0; j <= No_Lines; j++) {   // total number: No_Lines + 1
		kmax = mult;
		if (j == No_Lines)                  // proper number of data  
			kmax = resto;                   // for the last line 
		for (int k = 0; k < kmax; k++) {
	        char nb[20];
		    fscanf(from, "%s", nb);         // scan for space-delimited string
		    num = atof(nb);                 // convert string to a double value NIJE DOBRO
			i = j * mult + k;               // position in the cVector
		    setData(i,num);                 // assign data to a cVector
		}
	}
	// close the file
	fclose(from);
	if (Name != "cVECTOR")
		isNameGiven = true;}

//--- copy constructor
cVector::cVector(const cVector& original)
{
    data = new Complex[len = original.getLen()];
    assert(data != NULL);
    for(int i=0; i<len; i++) 
	  data[i] = original.data[i];
}

//--- destructor
cVector::~cVector()
{
	delete [] data;
}

//=============================== operator functions
//
//--- assignement operator
cVector& cVector::operator = (const cVector &original)
{
  if(this != &original) {
	delete [] data;
	data = new Complex[len=original.getLen()];
	assert(data != NULL);
	for(int i=0; i<len; i++) 
		data[i] = original.data[i];
	 Name = original.Name;
	 if (Name == "cVECTOR")
		 isNameGiven = false;
	 else 
		 isNameGiven = true;
  }
  return *this;
}
//--- assignement operator (set)
Complex& cVector::operator [] (int i)
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<len);
#endif
	return data[i];
}

//--- assignement operator (get)
const Complex& cVector::operator [] (int i) const
{
#ifdef ARRAYS_DEBUG
        assert(i>=0 && i<len);
#endif
        return data[i];
}

const Complex cVector::getData(int i) const
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<len);
#endif
	return data[i];
}

void cVector::setData(int i, Complex d)
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<len);
#endif
	this->data[i] = d;
	return;
}

//--- cVector B = self + v
cVector cVector::operator + (const cVector& v)
{
#ifdef ARRAYS_DEBUG
	assert(len == v.getLen());
#endif
        Complex zero(0.0,0.0);
	cVector sum(len, zero);
	for(int i=0; i<len; i++) 
		sum.data[i] = data[i] + v.data[i];
	return (sum);
}
//--- cVector B = self - v
cVector cVector::operator - (const cVector& v)
{
#ifdef ARRAYS_DEBUG
	assert(len == v.getLen());
#endif
        Complex zero(0.0,0.0);
	cVector sum(len, zero);
	for(int i=0; i<len; i++) 
		sum.data[i] = data[i] - v.data[i];
	return (sum);
}
//--- true if self == v
bool cVector::operator == (const cVector& v)
{
	bool test = false;
	if(len != v.getLen()) 
		return test;
	for(int i=0; i<len; i++) 
		if(data[i] != v.data[i]) 
			return test;
	// cVectors are equal
	return (test=true);
}

void cVector::Resize(int n)
{
	delete[]data;
	data = new Complex[len=n];
	assert(data != NULL);
}

void cVector::Zero(void)
{
	Complex num (0.0, 0.0);
	for(int i=0; i<len; i++) 
		data[i] = num;
	return;
}

const Complex cVector::Max(void) const
{
	Complex max = data[0];
	double norm = abs(max);   // modulus of a complex number = sqrt(x*x + y*y)
	for(int i=1; i<len; i++) 
		if(abs(data[i]) > norm) 
			max = data[i];
	return max;
}

const Complex cVector::Min(void) const
{
	Complex min = data[0];
	double norm = abs(min);
	for(int i=1; i<len; i++) 
		if(abs(data[i]) < norm) 
			min = data[i];
	return min;
}

const Complex cVector::Mean() const
{
	Complex sum (0.0, 0.0);
	for(int i=0; i<len; i++) 
		sum += data[i];
	return sum * (1./len);
}

Complex cVector::Dot(cVector& v)
{
	Complex s(0.0, 0.0);
#ifdef ARRAYS_DEBUG
       assert(len == v.getLen());
#endif
       for(int i=0; i<len; i++) 
	  s += data[i] * v[i];
       return s;
}

const Complex cVector::Modulus(void) const
{
	Complex s(0.0, 0.0);
	for(int i=0; i<len; i++) 
		s += data[i] * data[i];
	return sqrt(s);
}

//----------------- console I/O
//
void  cVector::rPrint()    // prints entries on a single line 
{
  cout << "cVector: " << Name << " = ";
  cout << "(";
  for(int i=0; i<len-1; i++) 
	  cout << data[i] << ",";
  cout << data[len-1] << ")" << endl;
  return;
}

void cVector::Print()
{
  cout << "cVector:  " << Name << endl;
  cout << "Length  = " << len  << endl;
  for(int i=0; i<len; i++) 
	  cout << data[i] << endl;
  return;
}

//------------------- File I/O
//
//--- read from predefined file
cVector cVector::Read()
{
	string FileName = Name + ".txt";
	cVector ans(FileName.c_str());
	return (ans);
}

cVector cVector::ReadFrom(const char* file)
{
	cVector ans(file);
	return (ans);
}

//--- write to a predefined file
void cVector::Write()
{
	int L = 20;   // default number of entries in a line
	if (len < L)
		L = len;
	string FileName = Name + ".txt";
	this->WriteTo(FileName.c_str(),L);
}

void cVector::WriteTo(const char* file, int L)
{
	int No_Lines;          // number of lines for ONE ROW 
	int resto;             // left-over in the last line for the row
	div_t Res;
	Res.quot = 1;  Res.rem = 0;
	if (len <= L)
		L = len;
	else if (len > L)             // if number of cols is > number of entries in a row
		Res = div(len, L);
	No_Lines = Res.quot;          // No_Lines with full number of mult entries
	resto    = Res.rem;           // ONE line with resto entries, all for just ONE ROW of data
	int J;
	Complex data;
	string VectorName = Name;
	//-- output to file 
	out.open(file,ios::out);
	assert(out != NULL);
	//--- write the HEADER BLOCK
	out << "!================== cVector(" << len << ") ====================" << endl;
	out << "N= " << VectorName << endl;
	out << "n= " << len <<  "  l= " << L  << endl;
	out << ":::" << endl;
	//--- write the DATA BLOCK
	int kmax;
	for (int j = 0; j <= No_Lines; j++) {
		kmax = L;
		if (j == No_Lines)
			kmax = resto;
			for (int k = 0; k < kmax; k++) {
				J = j * L + k;
				data = this->getData(J);
				out << data << "   ";
				out << this->getData(J) << "   ";
		}
		out << endl;
	}
	out.close();
}

//------------------ global cVector functions
//
Complex Dot(const cVector& v1, const cVector& v2)
{
  Complex s(0.0, 0.0);
#ifdef ARRAYS_DEBUG
  assert(v1.getLen() == v2.getLen());
#endif
  for(int i=0; i<v1.getLen(); i++) 
	  s += v1[i] * v2[i];
  return s;
}

cVector operator * (Complex x, const cVector& v)
{
        Complex zero(0.0,0.0);
	cVector ans(v.getLen(),zero);
	for(int i=0; i<v.getLen(); i++) 
		ans.setData(i,x*v[i]);
	return (ans);
}

cVector operator * (const cVector& v, double x)
{
        Complex zero(0.0,0.0);
	cVector ans(v.getLen(), zero);
	for(int i=0; i<v.getLen(); i++) 
		ans.setData(i,x*v[i]);
	return (ans);
}

Complex Modulus(const cVector& v)
{
	Complex s (0.0, 0.0);
	for(int i=0; i< v.getLen(); i++) 
		s += v.getData(i) * v.getData(i);
	return sqrt(s);
}

std::ostream& operator << (ostream& s, const cVector& v)
{
	s << "(";
	for(int i=0; i<v.getLen()-1; i++) 
		s << v.getData(i) << ", ";
	s << v.getData(v.getLen()-1)<<")"<<endl;
	return s;
}

/*==============================================
          cMatrix class definitions 
================================================*/

//--- constructors
//
cMatrix::cMatrix()
{
	m = new cVector[M_DEFAULT];
	assert(m != NULL);
    Complex zero(0.0,0.0);
	rows = cols = M_DEFAULT;                // cMatrix (2 x 2) with all zeroes
	cVector v(cols,zero);                        // cVector with all zeroes
	for(int i=0; i<rows; i++) 
		m[i] = v;
	Name = "cMATRIX";
	isNameGiven = false;
}

cMatrix::cMatrix(int r, int c, Complex d)
{
	m = new cVector[r];
	assert(m != NULL);
	rows = r;
	cols = c;
	cVector v(cols,d);
	for(int i=0; i<r; i++)         // cMatrix (r x c) with all elements set do d		
		m[i] = v;
	Name = "cMATRIX";
	isNameGiven = false;
}

cMatrix::cMatrix(int r, int c, string name)
{
	m = new cVector[r];
	assert(m != NULL);
	rows = r;
	cols = c;
	Complex d (0.0, 0.0);
	cVector v(cols,d);
	for(int i=0; i<r; i++)         // cMatrix (r x c) with all elements set do d		
		m[i] = v;
	Name = name;
	if (Name == "") {
		Name = "cMATRIX";
	    isNameGiven = false;
	}
	isNameGiven = true;
}

/**
** \brief Construct a cMatrix from a file
**
**    Format of the matrix file - it contains two parts: Header block and Data block
**    HEADER BLOCK:
**
**    (1) One or more optional COMMENT lines, where each comment line begins with '!'
**    
**    (2) Optional line with a NAME of the Matrix. If given, it must be written as:
**    N= (string) 
**
**    where the string in the name must begin with a letter (a-w,A-W).
**    If the Name is not given, default name MATRIX is assigned.
**
**    (3) Number of ROWS 'r', number of COLUMNS 'c' and number of elements in one line 'l'
**    r= (int)  c= (int)  l= (int)
**
**    (for small number of columns, c = l, i.e. each row is in one line)
**    NOTE - all three entries, r, c, and l must be given!!
**
**    The end of the Header part is given as 3 colons:
**    :::
**
**    DATA BLOCK:
**    numeric entris are given row by row (number of entries in a row is <= l)
**    (entries are separated by a blanko)
**
**    NOTE: For larger matrices (say, cols > 15-20), multiline input for each row is neccessary 
**    
**/

cMatrix::cMatrix(const char* file)
{
	const int BSIZE = 120;
    FILE *from;
    int row = 0, col = 0;       // number of rows and columns
    int mul = 0;                // number of entries in one line
 
	if ((from = fopen(file,"r")) == 0) {
		cerr << "\nCannot open " << file << endl;
		exit(1);
	}
    char buf[BSIZE], *cp, *cp2;
    int rfound = 0, cfound = 0, mfound = 0, nfound = 0, colonsfound = 0;  // control flags

    //---------------------------------------------------- read the HEADER BLOCK
    while( fgets(buf, BSIZE, from) ) {
		//--- read the comment lines
      if( ( cp = strpbrk(buf,"!")) != NULL )            
        *cp = '\0';                                     
		//---- read the Name of a Vector 
      if(( cp = strpbrk(buf,"N")) != NULL )          
        if (( cp2 = strpbrk(cp, "=")) != NULL )
          if (( cp = strpbrk( cp2, cAlphabet)) != NULL ) {
			  Name = cp;
			  isNameGiven = true;
			  nfound++;
            }   
		  //--- read the number of rows
	 if(( cp = strpbrk(buf,"r")) != NULL )
        if (( cp2 = strpbrk(cp, "=")) != NULL )
          if (( cp = strpbrk( cp2, "0123456789")) != NULL ) {
              row = atoi(cp);                           
              rfound++;                                 
            }
		  //--- read the number of columns
      if( ( cp = strpbrk(buf,"c") ) != NULL )
         if ( ( cp2 = strpbrk(cp, "=")) != NULL )
           if ( ( cp = strpbrk(cp2, "0123456789")) != NULL ) {
               col = atoi(cp);                          
               cfound++;                                
             }
		   //--- read the number of entries in a line
       if( ( cp = strpbrk(buf,"l") ) != NULL )
         if ( ( cp2 = strpbrk(cp, "=")) != NULL )
           if ( ( cp = strpbrk(cp2, "0123456789")) != NULL ) {
               mul = atoi(cp);                          
               mfound++;                               
             }
		   //--- termination of the Header block
     if ( strstr(buf,":::") != NULL ) {                 
          colonsfound++;
          break; //... out of "while" loop
        }
    }
    if ( !rfound || !cfound || !mfound || !colonsfound ) { // if the header block is finished
		cerr << "\nThe input file is NOT OK \n" << endl;
        exit(1);
    }
	if (!nfound) {                // no Matrix name is given
		Name = "cMATRIX";         // assign the default name
		isNameGiven = false;
	}
	//------------------------------------------------- make space for the cMatrix
	m = new cVector[row];              // array of row cVectors
	assert(m != NULL);
        Complex zero(0.0,0.0);
	cVector* v = new cVector(col,zero);     // each cVector with col elements
	assert(v != NULL);
	rows = row;
	cols = col;
	int I = 0;
	Complex data;
	// number of input lines for each cMatrix row (No_Lines) and 
	// number of remaining elements in the last line for a row (resto)
	// --- for one row of the cMatrix
	int No_Lines, resto;
	div_t Res;
	Res.quot = 1;  Res.rem = 0;
	if (col > mul)                // if number of cols is > number of entries in a row
		Res = div(col, mul);
	No_Lines = Res.quot;          // No_Lines with full number of mult entries
	resto    = Res.rem;           // ONE line with resto entries, all for just ONE ROW of data
    //-------------------------------------- read the DATA BLOCK
	int kmax;
	for (int i = 0; i < rows; i++) {
	    for (int j = 0; j <= No_Lines; j++) {
			kmax = mul;
			if (j == No_Lines)                 // proper number of data
				kmax = resto;                    // for the last line in a row  
		    for (int k = 0; k < kmax; k++) {     // read a row from a file
			    char nb[20];
			    fscanf(from, "%s", nb);          // scan for space-delimited string
			    data = atof(nb);                 // convert string to an COMPLEX value
			    I = j * mul + k;
			    v->setData(I,data);              // assign data to each row
		    }
		}  // end of data for one row
		m[i] = *v;
	}
	fclose(from);
	if (Name != "MATRIX")
		isNameGiven = true;
}

//--- copy constructor
//
cMatrix::cMatrix(const cMatrix& M)   
{
	int i;
	rows = M.getRows();
	m = new cVector[rows];
	assert(m != NULL);
	cols = M.getCols();
        Complex zero(0.0,0.0);
	cVector v(cols,zero);
	for(i=0; i<rows; i++) {
		v = M.getRowData(i);
		this->m[i] = v;
	}
}

//--- destructor
cMatrix::~cMatrix()
{
	delete [] m;
}

//--- assignement operator
cMatrix& cMatrix::operator = (const cMatrix& M)
{
	if(this != &M) {
		delete [] m;
		rows = M.getRows();
		cols = M.getCols();
                Complex zero(0.0,0.0);
		cVector v(cols,zero);
		m = new cVector[rows];
	        assert(m != NULL);
		for(int i=0; i<rows; i++) {
			v = M.getRowData(i);
			this->m[i] = v;
		}
	}
	return *this;
}

const Complex cMatrix::getElement(int i, int j) const
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<rows);
	assert(j>=0 && j<cols);
#endif	
	Complex val;
	val = this->m[i].getData(j);
	return (val);
}

void cMatrix::setElement(int i, int j, const Complex val)
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<rows);
	assert(j>=0 && j<cols);
#endif
	m[i].setData(j,val);
	return;
}

const cVector cMatrix::getRowData(int i) const
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<rows);
#endif
        return m[i];
}

const cVector cMatrix::getColData(int j) const
{
#ifdef ARRAYS_DEBUG
	assert(j>=0 && j<cols);
#endif
	int i, Rows = this->getRows();
	Complex val;
        Complex zero(0.0,0.0);
	cVector v(Rows,zero);
	for (i = 0; i < Rows; i++) {
		val = m[i].getData(j);
		v.setData(i,val);
	}
	return (v);
}

void cMatrix::setColData(int j, const cVector& v)
{
#ifdef ARRAYS_DEBUG
	assert(j>=0 && j<cols);
	int Vrows = v.getLen();
	assert(rows == Vrows);
#endif
	int i, Rows = this->getRows();
	Complex val;
	for (i = 0; i < Rows; i++) {
		val = v.getData(i);
		this->m[i].setData(j,val);
	}
}

void cMatrix::setRowData(int i, const cVector& v)
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<rows);
	int Cols = v.getLen();
	assert(cols == Cols);
#endif
	this->m[i] = v;
}

//--- Product of this cMatrix and a cVector
//--- cVector b = self * cVector v
cVector cMatrix::operator * (const cVector& v)
{
#ifdef ARRAYS_DEBUG
	assert(cols == v.getLen());
#endif
	cVector ans(rows); 
	Complex prod;
	for(int i=0; i<rows; i++) {
		prod = Dot(this->m[i],v);
		ans.setData(i,prod);
	}
	return (ans);
}

//--- Product of this cMatrix and a Scalar
//--- cMatrix B = self * Complex s
//
cMatrix cMatrix::operator * (const Complex s)
{
	cMatrix B(*this);      // resultant cMatrix: B = self * s
	Complex a_ij;
	for(int i=0; i<rows; i++) {
		for(int j=0; j<cols; j++) {
			a_ij = B.getElement(i,j);
			a_ij *= s;
			B.setElement(i,j,a_ij);
		}
	}
	return (B);
}

//--- Division of this cMatrix and a Scalar
//--- cMatrix B = self / double s
//
//cMatrix& cMatrix::operator / (Complex s)    
//{
//	if (s <= cTINY) 
//		exit(2);
//	cMatrix* B = new cMatrix(*this);      // resultant cMatrix: B = self / s
//	if (!B)
//		exit(2);	
//	Complex a_ij;
//	for(int i=0; i<rows; i++) {
//		for(int j=0; j<cols; j++) {
//			a_ij = B->getElement(i,j);
//			a_ij /= s;
//			B->setElement(i,j,a_ij);
//		}
//	}
//	return (*B);
//}

//--- Product of this cMatrix and cMatrix A
//--- cMatrix B = self * cMatrix A
//
cMatrix cMatrix::operator * (const cMatrix& A)
{
#ifdef ARRAYS_DEBUG
	assert(cols == A.getRows());
#endif
	Complex res;
	int Cols = A.getCols();
	cMatrix B(rows,Cols);  // resultant cMatrix: B = self * A
        Complex zero(0.0,0.0);
	cVector* v = new cVector(Cols,zero);       // columns of cMatrix A
	assert(v != NULL);
	// cMatrix product
	for(int i=0; i<rows; i++) {
		for(int j=0; j<Cols; j++) {
			*v  = A.getColData(j);     
			res = Dot(this->m[i],*v);    // scalar product
			B.setElement(i,j,res);
		}
	}
        delete v;
	return (B);
}

cMatrix cMatrix::operator + (const cMatrix& A)
{
	int i;
#ifdef ARRAYS_DEBUG
	assert(rows == A.getRows() && cols == A.getCols());
#endif
        Complex zero(0.0,0.0);
	cMatrix ans(A.getRows(),A.getCols(),zero);
	cVector* v   = new cVector(A.getCols(),zero);
	assert(v != NULL);
	for(i=0; i<A.getRows(); i++) {
		*v = (this->m[i] + A.getRowData(i));
		ans.setRowData(i,*v);
	}
        delete v;
	return (ans);
}

cMatrix cMatrix::operator - (const cMatrix& A)
{
#ifdef ARRAYS_DEBUG
	assert(rows == A.getRows() && cols == A.getCols());
#endif
        Complex zero(0.0,0.0);
	cMatrix ans(A.getRows(),A.getCols(),zero);
	cVector* v   = new cVector(A.getCols(),zero);
	assert(v != NULL);
	for(int i=0; i<rows; i++) {
		*v = (this->m[i] - A.getRowData(i));
		ans.setRowData(i,*v);
	}
        delete v;
	return (ans);
}

Complex cMatrix::RowSum(int i)
{
#ifdef ARRAYS_DEBUG
	assert(i>=0 && i<rows);
#endif
	Complex sum (0.0, 0.0);
	for(int k=0; k<cols; k++) 
		sum += m[i][k];
	return sum;
}

Complex cMatrix::ColSum(int j)
{
#ifdef ARRAYS_DEBUG
	assert(j>=0 && j<cols);
#endif
	Complex sum (0.0, 0.0);
	for(int k=0; k<rows; k++) 
		sum += m[k][j];
	return sum;
}

void cMatrix::Zero(void)
{
	Complex num (0.0, 0.0);
	for(int i=0; i<rows; i++) {
		for(int j=0; j<cols; j++)
			m[i][j] = num;
	}
	return;
}

void cMatrix::Identity(void)
{
#ifdef ARRAYS_DEBUG
	assert(rows == cols);
#endif
	Zero();
	Complex num (1.0, 0.0);
	for(int i=0; i<rows; i++) 
		m[i][i] = num;
	return;
}

void cMatrix::Diagonal(const Complex d)
{
#ifdef ARRAYS_DEBUG
	assert(rows == cols);
#endif
	Zero();
	for(int i=0; i<rows; i++) 
		m[i][i] = d;
	return;
}

void cMatrix::Diagonal(const cVector& v)
{
#ifdef ARRAYS_DEBUG
	assert(rows == cols);
#endif
	Complex d_i;
	Zero();
	for(int i=0; i<rows; i++) {
		d_i = v.getData(i);
		m[i][i] = d_i;
	}
	return;
}
/*------------- OUT
//===================== LU decomposition of cMatrix SELF
// (see GLOBAL ludcmp for details)
//  according to "NR in C"
//-----------------------------------------------------
//
cMatrix& cMatrix::calcLU (iVector& indx)
{
	cMatrix* ans = new cMatrix(*this);   // copy of SELF 
	if (!ans)                         
		exit(2);                      
	Complex d;
	int N = this->getRows();
	if (N != this->getCols())
		exit(2);
	//--- call ludcmp with original cMatrix ans
	ludcmp(*ans, indx, d);
	// LU decomposition is now cMatrix ans 
	return (*ans);
}

//=============================== solution of linear equations
//
//   equations: SELF * x = b  => solution x
//------------------------------------------------------------
//
cVector& cMatrix::solveEQ(iVector& indx, cVector& b)
{
	cMatrix* tmp = new cMatrix(this->getRows(),this->getCols());   
	if (!tmp)                      
		exit(2);                  
	cVector* ans = new cVector(b);       // copy of input cVector b
	if (!ans)                          // ans will hold the result
		exit(2);                       // after backsubstitution
//	double d;
	// control of input
	int n1,n2,n3;
	n1 = tmp->getRows();
	n2 = tmp->getCols();
	n3 = ans->getLen();
	if (n1 != n2 && n1 != n3) 
		exit(2);
	// call LU and solve
	*tmp = this->calcLU(indx);
	lubksb(*tmp, indx, *ans);
    // give solution
	return (*ans);
}

//============================ Inverse of a SELF cMatrix 
//
cMatrix& cMatrix::Inv()
{
	cMatrix* tmp = new cMatrix(*this);       // copy of SELF 
	if (!tmp)                          
		exit(2);                       
	int N = this->getRows();
	if (N != this->getCols())              // check if SELF is square
		exit(2);
	cMatrix* ans = new cMatrix(N,N);         // cMatrix will be Inverse of SELF
	if (!ans)                      
		exit(2);                   
	cVector* col = new cVector(N);   
	if (!col)                      
		exit(2);
	iVector* indx = new iVector(N);
	if (!indx)
		exit(2);
	Complex d;
//	algorithm from NR
	ludcmp(*tmp, *indx, d);
	for (int j = 0; j < N; j++) {
		col->setData(j,1.0);
		lubksb(*tmp, *indx, *col);
		for (int i = 0; i < N; i++) 
			ans->setElement(i,j,col->getData(i));
	}
	return (*ans);
}

//============================ Inverse of a SELF cMatrix 
//
cMatrix& cMatrix::operator ! ()
{
	cMatrix* tmp = new cMatrix(*this);       // copy of SELF 
	if (!tmp)                          
		exit(2);                       
	int N = this->getRows();
	if (N != this->getCols())              // check if SELF is square
		exit(2);
	cMatrix* ans = new cMatrix(N,N);         // cMatrix will be Inverse of SELF
	if (!ans)                      
		exit(2);                   
	cVector* col = new cVector(N);   
	if (!col)                      
		exit(2);
	iVector* indx = new iVector(N);
	if (!indx)
		exit(2);
	Complex d;
//	algorithm from NR
	ludcmp(*tmp, *indx, d);
	for (int j = 0; j < N; j++) {
		col->setData(j,1.0);
		lubksb(*tmp, *indx, *col);
		for (int i = 0; i < N; i++) 
			ans->setElement(i,j,col->getData(i));
	}
	return (*ans);
}
---------------- end OUT */
//============================== Transpose of a SELF cMatrix
//
cMatrix cMatrix::Tra()
{
	int N = this->getRows();
	int M = this->getCols();            // SELF may be N x M cMatrix
	cMatrix ans(M,N);     // ans will be transpose of SELF

	cVector* tmp = new cVector(M);       // row of SELF cMatrix
	assert(tmp != NULL);
	for (int i = 0; i < N; i++) {
		*tmp = this->getRowData(i);
		ans.setColData(i,*tmp);
	}
        delete tmp;
        // give solution
	return (ans);
}

//============================== Transpose of a SELF cMatrix
//
cMatrix cMatrix::operator ~ ()
{
	int N = this->getRows();
	int M = this->getCols();            // SELF may be N x M cMatrix
	cMatrix ans(M,N);      // ans will be transpose of SELF

	cVector* tmp = new cVector(M);       // row of SELF cMatrix    
	assert(tmp != NULL);
	for (int i = 0; i < N; i++) {
		*tmp = this->getRowData(i);
		ans.setColData(i,*tmp);
	}
        delete tmp;
        // give solution
	return (ans);
}
/*---------------- OUT
//=============================== Determinant of SELF cMatrix
//
Complex cMatrix::Det()
{
	int N = this->getRows();
	if (N != this->getCols())
		exit(2);
	iVector* indx = new iVector(N);
	if (!indx)
		exit(2);
	cMatrix* tmp = new cMatrix(N,N);   
	if (!tmp)                      
		exit(2);                  
	Complex D = 0.0;
	// call LU and solve
	*tmp = this->calcLU(*indx);
	for (int i = 0; i < N; i++)
		D += tmp->getElement(i,i);
	return D;
}

//=================== Division of SELF cMatrix with a cVector
//    i.e. Product SELF^(-1) * cVector v
//
cVector& cMatrix::operator / (cVector& v)
{
	cMatrix* tmp = new cMatrix(*this);  // copy of SELF 
	if (!tmp)                          
		exit(2);                       
	int N = this->getRows();
	if (N != this->getCols())         // check if SELF is square
		exit(2);
	if (N != v.getLen())              // check compatibility
		exit(2);
	cVector* ans = new cVector(N);       // result: SELF^(-1) * v    
	if (!ans)                          
		exit(2);
	*ans = !tmp * v;
	return (*ans);
}

//=================== Division of SELF cMatrix with a cMatrix
//    i.e. Product SELF^(-1) * cMatrix A
//
cMatrix& cMatrix::operator / (cMatrix& A)
{
	cMatrix* tmp = new cMatrix(*this);  // copy of SELF 
	if (!tmp)                          
		exit(2);                       
	int N = this->getRows();
	if (N != this->getCols())         // check if SELF is square
		exit(2);
	if (N != A.getRows())             // check compatibility
		exit(2);
	int M = A.getCols();
	cMatrix* ans = new cMatrix(N,M);    // result: SELF^(-1) * A    
	if (!ans)                          
		exit(2);
	//--- calculate
	*ans = !tmp * A;
	return (*ans);
}
------------- end OUT */
//----------------- console I/O
//
void cMatrix::Print() // print each row separately
{
	cout << "cMatrix Name = " << Name << endl;
	cout << "No. of Rows  = " << rows << " and Cols = " << cols << endl;
	for(int i=0; i<rows; i++) {
		 m[i].rPrint();
	}
}
void cMatrix::pPrint()
{
	for (int i=0; i<rows; i++) {
		for(int j=0; j<cols; j++) 
			cout << m[i][j] << " ";
		cout << endl;
	}
	return;
}

//------------------ file I/O
//
//!  read a Matrix from default file
cMatrix cMatrix::Read()
{
	string FileName = Name + ".txt";
	cMatrix ans(FileName.c_str());
	return (ans);
}

//!  read a cMatrix from a file
cMatrix cMatrix::ReadFrom(const char* file)
{
	cMatrix ans(file);
	return (ans);
}

//!  write Matrix to default file
void cMatrix::Write()
{
	int L = 20;                 // default number of entries in a line
	if (this->getCols() < L)
		L = this->getCols();
	string FileName = Name + ".txt";
	this->WriteTo(FileName.c_str(), L);
}

//!  write a cMatrix to a file
//
void cMatrix::WriteTo(const char* file, int L)
{
	int No_Lines;          // number of lines for ONE ROW 
	int resto;             // left-over in the last line for the row
	div_t Res;
	Res.quot = 1;  Res.rem = 0;
	if (cols <= L)
		L = cols;
	else if (cols > L)                   // if number of cols is > number of entries in a row
		Res = div(cols, L);
	No_Lines = Res.quot;                 // No_Lines with full number of mult entries
	resto    = Res.rem;                  // ONE line with resto entries
	                                     // all for just ONE ROW of data
	int J;
	Complex data;
	string MatrixName = Name;
	//-- output to file 
	out.open(file,ios::out);
	assert(out != NULL);
	//--- write the HEADER BLOCK
	out << "!================== cMatrix(" << rows << ", " << cols << ") ====================" << endl;
	out << "N= " << MatrixName << endl;
	out << "r= " << rows << "  c= " << cols << "  l= " << L  << endl;
	out << ":::" << endl;
	//--- write the DATA BLOCK
	int kmax;
	for (int i = 0; i < rows; i ++) {
		for (int j = 0; j <= No_Lines; j++) {
			kmax = L;
			if (j == No_Lines)
				kmax = resto;
			for (int k = 0; k < kmax; k++) {
				J = j * L + k;
				data = this->getElement(i,J);
				out << data << "   ";
			}
			out << endl;
		}
	}
	out.close();
}

//======================================= global cMatrix functions
//
//--- Multiply a cMatrix with a complex number
//
cMatrix operator * (const Complex x, const cMatrix& M)
{
	cMatrix ans(M.getRows(),M.getCols());

	for(int i=0; i<ans.getRows(); i++) {
		ans.setRowData(i, x * M.getRowData(i));
	  }
	return (ans);
}

cMatrix operator * (const cMatrix& M, const Complex x)
{
	cMatrix ans(M.getRows(),M.getCols());

	for(int i=0; i<ans.getRows(); i++) {
		ans.setRowData(i,x * M.getRowData(i));
	  }
	return (ans);
}
//--- Stream I/O
//
std::ostream& operator<<(ostream& s, const cMatrix& M)
{
	for(int i=0; i<M.getRows(); i++) 
		s << M.m[i];
	return s;
}
/*------------- OUT
//================================================= LU Decomposition (Crout's algorithm)
//
// A - input cMatrix of order nxn (destroyed afrer LU: A = LU) => output cMatrix
// indx - int cVector of order n => output icVector recording row permutations
//        (indx is neccessary for backsubstitution phase)
// d - output double (+/-1.0) depending if the number of row interchanges is even or odd
//        (d is not really neccessary)
// ludcmp is combined with lubksb to solve linear equations of invert a cMatrix
//
//               according to "Numerical Recipes in C"
//--------------------------------------------------------------------------------------
//
void ludcmp(cMatrix& A, iVector& indx, Complex& d)  
{
//	const double TINY=1.0e-20;
	int i,imax,j,k;
	Complex  sum;
	double big, dum, temp;
	int n = A.getRows();
	cVector vv(n);
	if (A.getRows() != A.getCols())
		exit(2);
	d = (1.0, 0.0);
	Complex a_ij, a_ik, a_kj, a_imaxk, a_jk, a_jj;
	for (i=0; i<n; i++) {
		big = 0.0;
		for (j=0; j<n; j++)  {
			a_ij = A.getElement(i,j);
			if ((temp = abs(a_ij)) > big) 
				big = temp;
		}
		if (big == 0.0) 
			cerr << "Singular cMatrix in routine LUDCMP" << endl;
		vv[i] = 1.0 / big;
	} // end for i = 0
	for (j=0; j<n; j++) {  
		for (i=0; i<j; i++) {
			a_ij = A.getElement(i,j);
			sum=a_ij;
			for (k=0; k<i; k++) {
				a_ik = A.getElement(i,k);
				a_kj = A.getElement(k,j);
				sum -= a_ik * a_kj;
			}
			A.setElement(i,j,sum);
		}
		big = 0.0;
		for (i=j; i<n; i++) {
			a_ij = A.getElement(i,j);
			sum = a_ij;
			for (k=0; k<j; k++) {
				a_ik = A.getElement(i,k);
				a_kj = A.getElement(k,j);
				sum -= a_ik * a_kj;
			}
			A.setElement(i,j, sum);
			if ((dum = vv[i] * abs(sum)) >= big) {
				big  = dum;
				imax = i;
			}
		}
		if (j != imax) {
			for (k=0; k<n; k++) {
				a_imaxk = A.getElement(imax,k);
				dum = a_imaxk;
				a_imaxk = A.getElement(j,k);
				a_jk = A.getElement(j,k);
				A.setElement(imax,k, a_jk);
				A.setElement(j,k, dum);
			}
			d = -(d);
			vv[imax] = vv[j];
		}
		indx[j] = imax;
		a_jj = A.getElement(j,j);
		if (a_jj == 0.0)
			A.setElement(j,j, cTINY);  
		if (j != (n-1)) {
			a_jj = A.getElement(j,j);
			dum = 1.0 / a_jj;
			for (i=j+1; i<n; i++) {
				a_ij = A.getElement(i,j);
				a_ij *= dum;
				A.setElement(i,j, a_ij);
			}
		}
	}
}

//============================================= BackSubstitution
//
//  cMatrix A is the LU decomposition of the system cMatrix
//     cMatrix A and icVector indx are output from ludcmp
//  cVector b is the right-hand side cVector
//  => the solution of A*x=b is contained in cVector b
//
//     cMatrix A and icVector indx are not modifed in lubksb
//     so, lubksb may be used for multiple right-hand sides
//--------------------------------------------------------------
//
void lubksb(cMatrix& A, iVector& indx, cVector& b)
{
	int i,ip,j;
	Complex sum;
	int n = A.getRows();
	Complex a_ij, a_ii;
	for (i=0; i<n; i++) {
		ip   = indx[i];
		sum  = b[ip];
		b[ip]= b[i];
		for (j=0; j<=i-1; j++) {
			a_ij = A.getElement(i,j);
			sum -= a_ij * b[j];
		}
		b[i] = sum;
	}
	for (i=n-1; i>=0; i--) {
		sum = b[i];
		for (j=i+1; j<n; j++) {
			a_ij = A.getElement(i,j);
			sum -= a_ij * b[j];
		}
		a_ii = A.getElement(i,i);
		b[i] = sum / a_ii;
	}
}

//============================== solution of linear equations
//
//   equations: A * x = b  => solution x
//-----------------------------------------------------------
//
cVector& solveEQ(cMatrix& A, cVector& b)
{
	cMatrix* tmp = new cMatrix(A);   // copy of the input cMatrix A 
	if (!tmp)                      // tmp will hold decomposition LU = A
		exit(2);                   // after LU: ludcmp
	cVector* ans = new cVector(b);   // copy of the input cVector b
	if (!ans)                      // ans will hold the result
		exit(2);                   // after backsubstitution: lubksb
	Complex d;
	// control of input
	int n1,n2,n3;
	n1 = A.getRows();
	n2 = A.getCols();
	n3 = b.getLen();
	if (n1 != n2 && n1 != n3) 
		exit(2);
	iVector indx(n1);
	// call LU and solve
	ludcmp(*tmp, indx, d);
	lubksb(*tmp, indx, *ans);
    // give solution
	return (*ans);
}

//============================ Inverse of a square cMatrix
//
//  A - input cMatrix (remains the same)
//   
//-------------------------------------------------------
//
cMatrix& Inv(cMatrix& A)
{
	cMatrix* tmp = new cMatrix(A);       // copy of input cMatrix A 
	if (!tmp)                          // tmp will hold decomposition LU = A
		exit(2);                       // to preserve the original cMatrix A
	int N = A.getRows();
	if (N != A.getCols())              // check if cMatrix A is square
		exit(2);
	cMatrix* ans = new cMatrix(N,N);     // cMatrix will be A_Inverse
	if (!ans)                      
		exit(2);                   
	cVector* col = new cVector(N);   
	if (!col)                      
		exit(2);
	iVector* indx = new iVector(N);
	if (!indx)
		exit(2);
	Complex d;
//	algorithm from "NR in C"
	ludcmp(*tmp, *indx, d);
	for (int j = 0; j < N; j++) {
		col->setData(j,1.0);
		lubksb(*tmp, *indx, *col);
		for (int i = 0; i < N; i++) 
			ans->setElement(i,j,col->getData(i));
	}
	return (*ans);
}
-------------------- end OUT */
//============================== Transpose of a cMatrix A
//
//  A - input cMatrix (remains the same)
//   
//-------------------------------------------------------
//
cMatrix Tra(const cMatrix& A)
{
	int N = A.getRows();
	int M = A.getCols();               // SELF may be N x M cMatrix
	cMatrix ans(M,N);     // ans will be transpose of SELF

	cVector* tmp = new cVector(M);       // row of SELF cMatrix    
        assert (tmp != NULL);
	for (int i = 0; i < N; i++) {
		*tmp = A.getRowData(i);
		ans.setColData(i,*tmp);
	}
        delete tmp;
    // give solution
	return (ans);
}
/*--------------- OUT
//=============================== Determinant of a cMatrix A
//
//  A - input cMatrix (remains the same)
//   
//---------------------------------------------------------
//
Complex Det(cMatrix& A)
{
	int N = A.getRows();
	if (N != A.getCols())
		exit(2);
	iVector* indx = new iVector(N);
	if (!indx)
		exit(2);
	cMatrix* tmp = new cMatrix(N,N);   
	if (!tmp)                      
		exit(2);                  
	Complex D = (0.0, 0.0);
	// call LU and solve
	*tmp = A.calcLU(*indx);
	for (int i = 0; i < N; i++)
		D += tmp->getElement(i,i);
	return D;
}
------------------- end OUT */

}   // end of namespace smv_libA

//============================================= end of file
